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Abstract 

A combination of ab initio , atomistic and finite element methods ( FEM ) were 
used to investigate the structures, energetics and lattice thermal conductance of 
grain boundaries for the ultra high temperature ceramic ZrB 2 . Atomic models of 
idealized boundaries were relaxed using density functional theory. Information about 
bonding across the interfaces was determined from the electron localization func- 
tion. The Kapitza conductance of larger scale versions of the boundary models were 
computed using non-equilibrium molecular dynamics. The interfacial thermal param- 
eters together with single crystal thermal conductivities were used as parameters in 
microstructural computations. FEM meshes were constructed on top of microstruc- 
tural images. From these computations, the effective thermal conductivity of the 
polycrystalline structure was determined. 
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I. INTRODUCTION 


Ultra high temperature ceramics (UHTC) are a class of materials with high melt- 
ing points, good mechanical properties and reasonable oxidation resistance. Among 
these materials, the metallic diborides, and especially ZrB 2 , have been the object of 
considerable study, both as pure materials and as constituents in composites. These 
materials are of interest for applications in extreme environments which require little 
or no oxidation or ablation. In particular, they are candidate materials for sharp 
leading edges and nosecaps of hypersonic aircraft as well as propulsion systems, re- 
fractory crucibles, among other applications [1-4]. Unlike most ceramics, UHTCs are 
distinguished by their high thermal conductivity. High thermal conductivity offers a 
number of advantages for high temperature applications. For example, high thermal 
conductivity can improve thermal shock resistance. It can also improve the efficiency 
of thermal radiation by rapidly redistributing thermal energy across available surfaces. 

The total thermal conductivity Ktot. of the diborides has significant contributions 
from both electronic and lattice carriers [5]. The electronic component K e can be 
estimated approximately from knowledge of the electrical conductivity a using the 
Weidemann-Franz (WF) empirical relation, n e = Kq uT, where k 0 is the Lorentz con- 
stant (2.45 x 10 _8 IU Wl ■ K ~ 2 ) and T is the temperature. The lattice contribution K p h, 
which results from phonon transport, cannot be measured directly and is usually in- 
ferred by subtracting n e from K tot . While not accessible to experimental measurement, 
the bulk lattice thermal conductivity as well as the lattice thermal resistance of grain 
boundaries can be determined directly from molecular dynamics (MD) simulations. 

For single crystal ZrB 2 , thermal conductivity measurements at room temperature 
have been reported as 140IU/ (m-K) in the basal direction and 100 W j (m-K) along the 
c-axis [6]. These single sample results did not include either a characterization of the 
defect distribution, which will reduce K to t , or specific estimates of K e and n p h- Poly- 
crystalline ZrB -2 has been more thoroughly studied. Previously, room temperature 
measurements have been reported as 55W/(m ■ K) for ZrB 2 where n e was estimated 
to be 33W/(m-K ) using the WF relation and K p h was given as 22 Wj (m-K) [7]. More 
recent results [8] quote values as high as 108 W/(m- K). The reduction of n tot relative 
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to single crystals comes primarily from the thermal resistance of grain boundaries. 
Variations of K to t with grain size can be substantial [7, 8] . While n e is expected to be 
relatively insensitive to grain boundaries due to the short electron mean free path, the 
lattice contribution K p h, which may be quite high in single crystals, can be reduced 
significantly due to grain boundary thermal resistance [7]. Different processing meth- 
ods can lead to very different grain sizes and grain boundary structures. Therefore, 
it is important to understand the effect of grain boundaries on the properties of these 
materials. 

In this paper, we use a multiscale approach to model structure-property relation- 
ships in polycrystalline ZrB 2 . We focus primarily on thermal properties; however our 
results also have implications for mechanical response. To start, ab initio methods 
were used to examine the structure and energetics of simple grain boundary struc- 
tures, namely coincidence and near coincidence tilt and twist boundaries. Relaxation 
of these interfaces leads to significant reconstruction of the atomic arrangement and 
also of the bonding across the interfaces. Thermal resistance of the boundaries was 
computed with nonequilibrium molecular dynamics using recently derived interatomic 
potentials for this material [9]. Previously, the lattice thermal conductivity of single 
crystal ZrB 2 was computed using these same potentials [10]. Thermal conductivity 
parameters from the atomistic simulations were used in microstructural calculations 
to determine the reduction in thermal conductivity due to the grain boundary net- 
work. Finite element method ( FEM ) computations were performed where FEM 
meshes were constructed directly on scanning electron microscopy ( SEA4 ) images. 
From these computations, the “effective thermal conductivity” of the polycrystalline 
microstructure could be determined. 

II. METHOD 

Ab initio computations were performed in the context of density functional theory 
(DFT) utilizing the functional of Perdew, Burke, and Ernzerhof (PBE) [11], This 
functional is known to give good results especially for solids. Ab initio methods were 
used to relax the grain boundary structures and to evaluate their energetics. Ab 
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initio results were compared to results obtained from the interatomic potentials to 
help assess the accuracy of those potentials. The plane wave code VASP was used 
for all calculations [12]. The Projected Augmented Wave (PAW) potentials were 
employed. All results were converged with respect to k-point sampling and plane 
wave energy cutoff E c . Because of the large unit cells considered, small h-space 
meshes were sufficient. Typically, 2x2x2 meshes were used for larger systems and 
6a;6a:6 meshes for smaller unit cells. Energy cutoffs were typically about 400eV. 

Information about bonding was obtained from the electron localization function 
(ELF) [13, 14]. The ELF gives the probability of finding an electron near a reference 
electron at a given point and with a fixed spin. In this way, it is useful for identifying 
covalent bonds. The ELF is defined in terms of a dimensionless localization ratio 
Xa( r ) which measures localization with respect to the uniform electron gas 


ELF(r) 


1 

1 + Xl(r)' 


( 1 ) 


The ELF takes values in the range 0 < ELF < 1 where ELF = 1 corresponds to 
perfect localization and ELF = 0.5 gives the electron gas. Xa( r ) is computed directly 
from DFT quantities. The ELF can be contrasted against the electron charge density 
p(r) which is a single electron quantity whereas the ELF gives information about the 
two-body distribution. 

Nonequilibrium molecular dynamics (MD) simulations were performed with the 
LAMM PS package [15] and using the interatomic potentials for ZrB 2 . The method 
of Muller-Plathe (MP) was used to compute interfacial thermal conductances [16 
18]. Periodic unit cells with a long direction normal to the boundaries were created. 
Periodicity results in two boundaries per cell. NVE simulations were performed 
in which a “hot” region and a “cold” region were created on opposite sides of the 
boundaries by exchanging atoms between the two regions: the atom with the greatest 
kinetic energy in the cold region is swapped with the lowest kinetic energy atom in 
the hot region. After a steady state is established, a temperature gradient between 
the regions, and across the interfaces, resulted where the heat flux Q is given by 


Q = 



( 2 ) 
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where the sum is over the number of exchanges, t is the total simulation time, A is the 
cross-sectional area normal to the direction of heat flow, m is the atomic mass and v % 
is the velocity of the exchanged atom. Near the sink/source regions, the temperature 
profile is nonlinear. For a homogeneous system however, i.e. without any interfaces, 
a linear region between the sink/source will develop whose slope obeys Fourier’s Law. 
In the presence of an interface, a grain boundary, for example, a discontinuity in the 
temperature profile will appear at the interface. The temperature discontinuity and 
the exchanged heat flux are related by the interfacial or Kapitza conductance <jk 

Q = a y<-AT. (3) 

where the inverse of the Kapitza conductance is the Kapitza resistance Rk = 1 /ctr- 

Simulations were performed using the recently derived ZrB 2 interatomic potentials 
of Daw, Lawson and Bauschlicher (DLB) [9]. Two sets of potentials were reported: 
“Pot 1” and “Pot 2”. Subsequently, in lattice thermal conductivity simulations [10], 
Pot 2 was found to give an normal conductivity that was greater than the in-plane 
value, a result that contradicted the experimental ordering. Pot 1 however gave 
generally good results including the correct ordering. For that reason, in this work, 
we focus on Pot 1 only. 

Simulations were performed as follows. Normal to the interfaces (^-direction), the 
unit cells had dimensions on the order of lOOnm to minimize reflections from the 
boundaries. Results were only weakly dependent on the lateral dimensions, however, 
which was typically (~ 1 — 5 nm). Convergence studies for both the normal and 
lateral dimensions were performed. Initial velocities were generated from a gaussian 
distribution. A NPT simulation was run at the given temperature for lOOps to equi- 
librate the system. Time steps of 1/s were typically used. Next, particle swapping 
according to the Muller-Plathe algorithm was performed; steady state was typically 
well-established after 500ps. Temperature profiles (T vs z) were determined by divid- 
ing the system into narrow slices perpendicular to the long direction of the system, 
i.e. parallel to the grain boundaries. Within each slice, the local temperature was 
evaluated. Steady state temperature profiles were determined by averaging the tem- 
peratures in each slice over runs of 4ns. The temperature discontinuity across the 
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interface was determined by fitting linear functions on opposite sides of the disconti- 
nuities and evaluating the temperature drop at the interface. 

Finite element computations were performed using the open source code OOF 2 [19] 
and the commercial code M SC Marc [20]. Finite element meshes were constructed 
directly on an SEM image of the microstructure of ZrB 2 [3]. OOF 2 contains image 
processing and mesh generation algorithms especially suited for imaged based FEM. 
MARC was used to perform the FEM computations because of its more extensive 
computational capabilities. FEM meshes were created for both the grains and the 
boundaries, each of which was treated as a separate material with distinct thermal 
conductivities. Thus, the grain boundaries are treated as an “interphase”. FEM 
interface computations are often performed using boundary elements however OOF 2 
does not currently have that capability. The two descriptions are expected to be 
equivalent. 

The effective thermal conductivity K e ff of the full microstructure (grains + bound- 
aries) was evaluated in three different ways. First, a transient thermal analysis was 
performed where a thermal flux was applied to one side of the system while the other 
sides were maintained as adiabatic boundaries. The temperature of the opposing side 
of the system (the “backface”) was determined as a function of time. The K e ff of 
the full polycrystal was evaluated by performing a second transient thermal anal- 
ysis on an equivalent, homogeneous, reference material. The thermal conductivity 
of the reference was tuned to give a backface temperature trace that matched the 
polycrystalline result. 

Next, we performed a steady state thermal analysis using two different sets of 
boundary conditions. First, a uniform temperature gradient ( UTG ) was applied by 
fixing the temperature of opposite sides of the model. A transient thermal analysis 
was performed until the system reached steady state. Second, we applied a uniform 
heat flux ( UHF ) to opposites sides of the system until a steady state was reached. 
For both cases, UTG and UHF , 

< q>= -Keff < VT > (4) 

where q is the heat flux, K e ff is the effective thermal conductivity, VT is the tem- 
perature gradient and the brackets represent volume averages over the system. Note 
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for UTG, we fix VT and measure < q > after steady state is achieved. For UHF, 
we do the converse. Both methods were performed with boundary conditions ap- 
plied in both the vertical and horizontal directions as well, ft has been shown that 
kuhf < k < k-utg [ 21 ]. 

A fundamental question in evaluating the effective properties of heterogeneous, 
microstructural models is whether the model is representative of the bulk system 
[22, 23]. It has been suggested that a model can be considered representative if its 
response to different boundary conditions, such as UTG and UHF, is the same. Thus, 
any microstructure above some minimal size will yield the same effective properties 
independent of boundary conditions, i.e. the bounds on k are very tight. If the 
model is too small, i.e. it is not representative, different boundary conditions will 
give different effective properties. 


III. STRUCTURAL MODELS 

The single crystal structure of ZrB 2 is the AIB 2 - type, designated as U32 with 
space group symmetry P6/mmm. This is a layered arrangement with alternating 
planes of pure Zr and pure B atoms. The Zr atoms are transition metals and are 
considerably larger than the B atoms both in terms of atomic radii and also mass. 
The Zr atom has a mass of 91.2 g/mol and an atomic radii of 1.6A while Boron is a 
relatively small atom with an atomic radius of 0.9A and an atomic mass of 10.8 g/mol. 
The Zr planes have a hexagonal close packed structure while the B planes are open, 
hexagonal with six membered rings resembling graphite. The layers are situated such 
that each Zr atom lies directly above and below the centers of 6-membered B rings 
in the adjacent planes. The primitive cell of these structures contains one formula 
unit (one Zr and two Bs) and the crystal has simple hexagonal symmetry (Deh). The 
lattice constant “a” gives the metal-metal atom distance within a plane and the “c” 
lattice constants give the metal-metal atom distance in alternate planes. The a and b 
axes are symmetry-equivalent, in-plane directions while the c axis is normal to planes. 

Structural models for ZrB 2 grain boundaries are considered for four situations: 
two twist boundaries and two tilt boundaries, one for each of the two independent 



axes, a and c. Hexagonal crystals with irrational ratios of their lattice parameters 
(c/a) 2 only have exact coincidence site lattice (CSL) boundaries corresponding to 
rotations about the c-axis [24-26]. We consider such CSL boundaries formed by 
twists and tilts about that axis. However, near coincidence boundaries about the 
a-axis can be formed by approximating (c/a) 2 as an rational number ( m/n ). This 
approximation will introduce a strain into the system. Smaller (n, m) values result 
in smaller, computationally efficient unit cells, but with larger strain. For larger 
values, the converse is true. For our structures, the strain was typically less than 1%. 
To accommodate ab initio computations, periodic cells were used, containing two 
boundaries per cell. With small cells, we expect nontrivial interaction between the 
boundaries. However we should be able to obtain reliable information about relative 
structural and energetic differences. 

For ceramics such as ZrB 2 , grain boundaries are often disordered and may contain 
many impurities. However, such boundaries are usually the result of particular pro- 
cessing methods. Cleaner and more crystalline boundaries are expected to have low 
interfacial energies and therefore their formation should be favored with improved pro- 
cessing methods. Further, such boundaries should have greater interfacial adhesion 
and reduced thermal resistance, and therefore are more desirable. Such boundaries 
are routinely observed in other ceramics such as alumina and zirconia [27, 28]. In this 
paper, we focus on these relatively simple boundaries both because of their compu- 
tational efficiency and also because of their expected superior properties. We expect 
grain boundary thermal properties to be somewhat insensitive to the exact atomic 
structure of the boundary. Therefore, it is not necessary to find the absolute, lowest 
energy boundary structure, which can be a highly nontrivial task as has been shown 
in recent work [29, 30] and is beyond the scope of this paper. 

A. C-tilt 

The first boundary we consider is a E7 symmetric tilt around the c-axis. This 
exact CSL boundary is shown in Fig 1. As a shorthand we refer to this boundary as 
“c-tilt”. Because the Boron sublattice in ZrB 2 is graphitic, we propose a c-tilt grain 
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boundary structure in analogy with graphite. Related structures have been recently 
proposed for graphene [31, 32], The boundary consists of a sequence of five and 
seven (5-7) membered Boron ring units that are separated by hexagonal rings. The 
distance between these units is related to the tilt angle where lower angles correspond 
to a larger separation between the 5-7 pairs. The smallest such structure is the £7 
structure. To complete the structure for ZrB 2 , Zr atoms are positioned over/above 
all Boron rings including the 5 and 7-membered rings. This structure can also be 
viewed as an array of edge dislocations with a horizontal Burger’s vector where the 
five membered rings represent the extra plane of atoms. 

B. C-twist 

The second boundary is a £7 symmetric twist about the c-axis. We call this 
boundary “c-twist” and it is displayed in Fig 2. This boundary is the simplest of 
the ones we will consider and shows the least reconstruction. The interface structure 
has a layer of Zr and a layer of B shifted relative to each other. Thus, across the 
interface, the Zr atoms are not centered directly above and below Boron 6-membered 
rings. This twist results in a small degree of crumpling of the two layers, mainly in the 
Boron layer. However, large scale reconstruction does not occur because the intralayer 
interactions (metallic Zr and covalent B) have not been disturbed significantly. The 
interfacial layers also are further stabilized by the compounded effect of additional 
layers away from the boundary. 

C. A-tilt 

The third structure we consider is a 90°, near coincidence, asymmetric tilt bound- 
ary about the a-axis, which we call “a-tilt”. This model contains two boundaries, 
which we designate Left and Right, as shown in Fig 3. This configuration has atomic 
layers oriented perpendicular to each other. With respect to the Boron layers, the 
Left and Right boundaries are mirror images of each other. Namely, flat, vertical, 
graphitic sheets of Borons (the outer sections of Fig 3) interact with Boron edge con- 
figurations (the inner section of Fig 3). This mirror symmetry is broken however by 
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Zr layers inserted into the two boundaries. A rotated view of this structure can be 
seen in Fig 7. 

For the Left boundary, the inserted Zr layers match the middle section Boron 
layer. However, a deficit of Zr atoms now exists to match the Boron layer to the 
left. Recall the Zr atoms want to be positioned inside a Boron 6-membered ring. 
This is only possible for the top and bottom Zr atoms in the Left boundary of Fig 7, 
leaving the middle two Zr atoms to share three Boron rings. As expected, the Boron 
rings with the perfect Zr match remain nearly flat, however the next, neighboring 
ring crumples slightly due to the two, off-center, middle Zr atoms. The next ring, 
on the other hand, is flat again as a result of a balancing between the two, off-center 
Zr atoms. This is somewhat surprising since there are not two Zr atoms centered on 
opposite sides of this ring. Thus, we find an alternating structure of flat and crumpled 
Boron rings in the Left interface. 

For the Right boundary, the Zr sheet matches perfectly the Boron layer to its 
right. However, there are too many Zr atoms in this layer to match the edge Boron 
layers on its left. These extra Zr atoms produce considerable crowding of the middle 
section Borons as can be seen in Fig 7. In fact, we see that every fourth Zr atom lines 
up “eye-to-eye” with a edge-on Boron sheet. These “edge-on” Zr atoms compress the 
other three Zr so that they cannot line-up with the Zr sheets in the middle section. 
This Zr compression has an effect on the Boron sheets which are bent up or down at 
their edges. 

D. A-twist 

The final boundary is a 90°, near coincidence, asymmetric twist about the a-axis, 
which we designate “a-twist”. There are important similarities between a-twist and 
a-tilt. Fig 4 shows the periodic model for a-twist with two boundaries, denoted Left 
and Right. Both Left and Right boundaries have two sets of 90° rotated edge-on 
Boron planes interacting with each other. 

With the Left boundary, the plane of Zr atoms matches the layering sequence of 
the middle section of the structure. However, the edge-on Boron planes to the right 
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of the interface want a Zr fitted into each well of the zig-zag Boron configuration. On 
the Left boundary, therefore, there is one Zr atom too few per repeat unit, leaving 
an effective Zr vacancy in the center of the boundary. This results in a complicated 
4-membered Boron ring configuration. Above and below this ring, Zr atoms are 
pushed away, resulting in the bending of both Zr and B planes near that boundary. 

The Right interface is different in that the plane of Zr atoms matches precisely 
the edge-on Boron planes to the left of the boundary. However, there is an extra 
Zr atom (actually a column of Zr going into the page of the figure) in this layer 
relative to the middle section where the Zr atoms want to be positioned between the 
Boron layers. The extra Zr atom, however, creates the same situation as a-tilt in 
that the Zr atom is effectively attached to the end of a Boron plane. Zr atoms above 
and below this atom are pushed away causing bending of the middle planes as they 
approach the boundary. We can see in the Boron plane that is edge-on to the extra 
Zr atom that the Boron rings are compressed along that edge. A rotated version of 
this exact situation is seen in the Left boundary where a line of edge-on Borons has 
been pushed back behind the main part of the boundary. This column of Borons has 
been compressed by an extra column of Zr that is vertical to the figure on the Left 
boundary. 

While a-twist and a-tilt have similar features, a-twist appears more disordered than 
a-tilt. This is a result of how the Boron planes interact with each other across their 
respective boundary. For a-tilt, the Boron planes on the left of the Left boundary and 
on the right of Right boundary are flat, chemically unreactive with no dangling bonds. 
The middle section of a-tilt has exposed edge-on Borons that could be reactive, but 
they are separated from the unreactive sheets by layers of large Zr atoms. For a-twist, 
the Boron sheets on both sides of both boundaries have the more reactive, edge-on 
configuration. This permits nontrivial Boron-Boron bonding across the interface and 
therefore a denser boundary. We will have more to say about this when we examine 
the electronic structure. 
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IV. ELECTRONIC STRUCTURE 


The electronic structure of single crystal ZrB 2 , and the related material HfB 2 , 
was discussed in detail in a previous publication [33]. Here we summarize some of 
those results. The electronic configuration for Zr is [Kr]f>s 2 Ad 2 . Thus, this species has 
two “s” and two “d” valence electrons to donate to the material. Boron’s electronic 
configuration is [He]2s 2 2p 1 , giving it one less valence electron than Carbon. This 
electron “deficiency” has a significant impact on the properties of pure Boron whose 
ground state structure is still unresolved. In pure materials, Boron is known to form 
strong covalent bonds. 

Despite its simple atomic structure, ZrB 2 displays all three major electronic bond- 
ing motifs. Namely, the bonding in the Zr planes is metallic; the bonding in the B 
planes is covalent; and the bonding between layers is ionic. In the metallic Zr layers, 
the electron localization function {ELF) is diffuse and relatively non-localized. ELF 
accumulation points however were found in the Zr planes positioned at the center 
of the Zr triangles. Thus these points formed six-member rings in the Zr plane. 
Interestingly, the accumulation points appear above/below Boron atomic sites in the 
neighboring planes. It is possible that the interaction between the Zr electron den- 
sity at the ELF accumulation points and the Boron density determines the nature of 
the interlayer bonding. In the B planes, the electrons are well-localized into covalent, 
highly directional bonds between neighboring B atoms. This is reflected in high ELF 
values between neighboring Boron atoms. Between layers, significant charge trans- 
fer was found as electron “rich” Zr layers donated electronic charge to the electron 
“deficient” B layers. 

The electronic structure of grain boundaries can provide important information 
that impact interface properties [34, 35]. The formation of an interface may result 
in significant reconstruction of the boundary atomic structure. The accompanying 
changes to the electronic structure can result in the formation of localized electronic 
states at the boundary as well as significant modifications to the bonding across the 
boundary. A range of scenarios is possible from the formation of dangling bonds to 
strong covalent bonding across the interface. The nature of the interfacial bonding 
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will affect boundary properties such as adhesion, formation energies, and thermal 
conductance, among others. 

A. C-tilt 

In Fig 5 we show the ELF for c-tilt. C-tilt is one of the least disrupted boundaries. 
The material is still layered in the same way as with the single crystal, but at the 
boundary, 5-7 ring pairs are inserted which results in the different crystal orientation 
on either side of the boundary. The ELF clearly shows the same bonding motifs 
survive the formation of this boundary, namely there is covalent bonding in the Boron 
planes including among the 5-7 rings. Region of high ELF are seen between all Boron 
pairs regardless of the ring type (5-,6-,7-membered rings). Interesting, the 5-7 rings 
remain flat without crumpling. Because of the strong covalent bonding across this 
interface, we expect especially strong interfacial adhesion and thermal conductance 
for this boundary. 

B. C-twist 

The bonding structure for c-twist is relatively simple. The crystal structure on 
either side of the boundary remains intact except at the interface where a Boron and 
a Zirconium layer have been shifted relative to each other. We display the ELF for 
the boundary in Fig 6. In particular, we consider the ELF in the plane on the Zr 
side of the boundary. This particular plane is sandwiched between two Boron planes. 
In the single crystal, the Boron planes are mirror images and the Zr ELF forms 
accumulation points in an open hexagonal lattice that corresponds to the locations of 
the Borons in the two neighboring planes. In the case of c-twist however the Boron 
in the two neighboring planes are shifted relative to each other. Therefore, the ELF 
accumulation points in the Zr plane no longer line up with the Boron atoms across 
the interface. The implication of this is not completely clear, but it is expected that 
this mismatch will weaken the ionic bonding across the interface. 
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C. A-tilt 


In Fig 7, the ELF for a-tilt is shown in a plane crossing the boundary where 
the orientation has been rotated with respect to Fig 3. Fig 7 shows a contour slice 
where blue is high ELF (> 0.5), green is low ELF (< 0.5) and white is close to 0.5. 
As with Fig 3, this cell has two asymmetric grain boundaries, Left and Right. The 
middle section of the Figure has horizontal layers while the left and right sections 
have vertical layers, reflecting the 90° tilt of the two boundaries. The ELF plane 
displayed in Fig 7 contains only Boron atoms on both sides of the interfaces. The Zr 
atoms in the Figure are in front of the contour plane. The plane that is displayed 
has the maximum ELF at the interface. High ELF in Z r is associated principally 
with the Boron sublattice where covalent bonding exists between Boron atoms. Thus, 
pockets of high ELF result from pairs of Boron atoms and indicate possible covalent 
bonding across the boundary. 

As discussed previously, the Left boundary is an interface between a flat, vertical 
Boron layer and the neighboring Zr layer to the right. The Zr layer in this interface 
has missing atoms with respect to the Boron layer, and therefore not enough Zr 
atoms match with the Boron rings to the left. The bonding across this boundary is 
expected to be ionic due to the Zr — B interaction However, we see evidence from the 
ELF plot of some non-trivial electron localization inside the interface involving pairs 
of B atoms. In Fig 7, the vertical Boron sheet to the left has three pockets of high 
ELF which line up with five pockets of high ELF from the horizontal sheets. These 
Boron ELF clouds are distorted slightly and may enhance the interaction across the 
interface. On the Right boundary, we have too many Zr atoms, and in fact, one 
column of the middle Zr (into the page of Fig 7) is partially surrounded by an “edge- 
on 1 ' Boron ring as opposed to hovering above or below. The Zr atoms above and 
below this edge-on Zr are pushed up and down causing corresponding bending of 
the Boron planes. From the ELF plot, pockets of localized electrons are also seen 
associated with B pairs atoms across the boundary. In this case, some of the B 
pairs have been distorted relative to their position on the Left boundary. The Right 
boundary is more complex than the Left boundary since there are Zr-Zr interactions 
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in addition to Zr-B and B-B interactions. 


D. A-twist 

In Fig 8, we show an ELF contour plot for a-twist; the plane displayed has the 
maximum ELF across the interfaces. As discussed previously, a-twist is an asym- 
metric configuration and therefore the Left and Right boundaries are different. Due 
to the reactive edge-on Boron plane interacting at both interfaces, a-twist has com- 
plex bonding across its boundaries. A central exotic feature of this structure is the 
horizontal Boron plane in the center of the middle section of Fig 8. This Boron plane 
intersects the two boundaries in very different ways. At the Left boundary where 
there are missing Zr atoms, 4-membered covalent Boron rings form. The ELF shows 
that there is strong covalent bonding in this ring that reaches across the Left interface. 
At the Right boundary, the same Boron plane encounters an extra Zr atom at that 
interface. As as result, this particular Boron plane is compressed relative the other 
Boron planes. There is strong covalent bonding across both Left and Right interfaces, 
resulting from 6-membered rings twisted at 90° relative to each other across the in- 
terfaces, but still covalently bonded to each other. In addition, most of the Zr atoms 
in both boundaries are surrounded by Boron structures on five sides, instead of the 
single crystal configuration with Borons only above and below the Zr. Thus, there is 
non-trivial and complex electronic structure between these Zr and their Boron rich 
caged structures. 

V. GRAIN BOUNDARY ENERGETICS 

The grain boundary interfacial excess energy was defined as 

7 = {Eqb — nE bu i k )/2A (5) 

where E GB is the total energy of the grain boundary cell, E buik is the total energy 
of a formula unit of single crystal bulk material, n is the number of formula units 
to match the number of atoms in the grain boundary cell, and A is the area of the 
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interface. The factor of “2” is because there are two boundaries per periodic cell. 
The volume of the interface defined was as 

Ax = ( zqb ~ Zh u ik)/2 (6) 

where zqb is the length of the relaxed cell in the direction normal to the boundary 
and Zb u ik is the length of an equivalent cell with no boundary. 

These quantities were calculated for unit cells corresponding to the structures 
shown in Fig 1-4 using ab initio computations and also using the ZrB^ interatomic 
potentials. The results are shown in Table I. With such small cells (from 72 atoms 
for a-tilt to 1344 atoms for a- twist), significant interactions is expected between the 
two boundaries in each cell. As previously discussed, c-tilt and c-twist are exact 
CSL boundaries whereas a-tilt and a-twist are near CSL structures. The near CSL 
interfaces have a small amount of artificial strain (less than 1%). To maintain these 
structures, their volume was fixed during the optimization and only the ionic positions 
were allowed to relax. Because of the fixed volume, A- will necessarily be zero for 
these cases. Both of these factors, interface interactions and unit cell strain, will 
affect our results in absolute terms however relative differences can be considered. In 
addition, comparison with ab initio results will help access the accuracy of these new 
potentials. 

In Table I, the relaxed and unrelaxed (in parenthesis) interfacial energies are given. 
The amount of relaxation in the energy is related to the degree of reconstruction of the 
structures. C-tilt and c-twist have the lowest interfacial energies with 153 meV/A 2 
and 157 meV/A 2 , respectively, compared to 227 meV/A 2 for a-tilt and 212 meV/A 2 
for a-twist. A-tilt and a-twist have considerably more reconstruction compared to 
c-tilt and c-twist as indicated by the significant relaxation in these structures. It is 
somewhat surprising that a-tilt has a higher energy than a-twist since a-twist is more 
disordered. However, the significant covalent bonding across the a-twist interfaces 
may be a factor in reducing its energy relative to a-tilt. The results for c-tilt and c- 
twist are consistent with the intuition that they are the least reconstructed structures 
and therefore should be the easiest to form. They are also low E CSL boundaries and 
therefore are expected to have low energies. An ab initio number for A, can only be 
evaluated for c-twist and that value is 0.29A. 
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Interfacial energies from the interatomic potential are very similar among the 
boundaries, with a variation of a few percent. Results are consistently lower for 
the potentials compared to the ab initio numbers, in the case of a-tilt and a-twist by 
a factor of 2X. The highest energy from the potentials is 118 rneV/A 2 for a-twist and 
the lowest is 107 rneV/A 2 a-tilt. Given the similarity of these numbers, it is difficult 
to establish a clear ordering. However, the degree of relaxation is consistent with the 
ab initio results with a much small amount of relaxation for c-tilt and c-twist com- 
pared to a-tilt and a-twist. The results for for c-twist agrees very well with the ab 
initio number. It is clear however that the DLB potential does better for unrelaxed 
energies than relaxed. In fact, the expected ordering of the relaxed boundaries, as 
suggested by the ab initio results, are not well reproduced by DLB. 

VI. THERMAL RESISTANCE 

In Table I, we also present values for the thermal conductance of each boundary. 
These calculations were all performed at 300K. As discussed previously, the values 
were evaluated by creating a thermal gradient across the boundaries and evaluating 
the magnitude of the discontinuity at the interfaces. Thermal profiles for the four 
boundaries are shown in Fig 9. Unit cells with long dimensions normal to the bound- 
ary were used. These cells were typically more than 200nm in length to minimized 
reflections. The longest cell had of length of 248nm for c-twist. A-twist which had an 
unusually large lateral cross-section was the smallest at 109nm in length. Note that 
c-tilt and c-twist are symmetric boundaries whereas a-tilt and a-twist are asymmet- 
ric. Therefore, if the two boundaries in a unit cell are symmetric, they will both have 
same a k- hi the asymmetric case, they can be different as indicated in the Table. 

Some degree of correlation is expected between the interfacial energy 7 and the 
thermal conductance a K . In particular, boundaries with low 7 should have high 
<7 k . The variation among the <jk values however is much greater than the variation 
among 7 values from the potentials. Thus, establishing a clear relation is not obvious. 
However, the ok values do correspond to our intuition for these structures. C-tilt 
has the highest conductance with 1.77 GW / (K -m 2 ). The fact that c-tilt should have 
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a high thermal conductance is not surprising since the structural disruption of that 
boundary is relatively small and the strong covalent bonding across the boundaries 
should permit a smooth flow of thermal energy. A-twist and c-twist have very similar 
conductances with 0.55/0.53 GW/(K ■ m 2 ) and 0.58 GW j{K ■ m 2 ). Their 7 values 
from the potentials are very similar however a much larger gap is seen in their ab initio 
7 values. A-tilt has the lowest conductances with values of 0.38/0.31 GW/(K ■ m 2 ) 
for the Left and Right boundaries. The a-tilt 7 value is similar to a-twist at the ab 
initio level. However, we expect a-twist to have a higher thermal conductance than 
a-tilt because of the significant covalent bonding across that interface. 

It is interesting that c-tilt which has a 7 values on par with c-twist at the ab initio 
level has a much larger value for the thermal conductance Ok ■ This may be due to 
the fact that thermal conductivity values are higher in-plane compared to normal 
directions. Therefore, in-plane interfacial conductances may be expected to be higher 
than for normal interfaces. 

VII. MICROSTRUCTURAL MODELING 

The effect of grain boundary networks on the bulk thermal conductivity of poly- 
crystal ZrB 2 is now considered. To begin, a Brick Layer Model (. BLM ) is used to 
obtain an estimate for the reduction of thermal conductivity due to interfacial ther- 
mal resistance [36]. In the BLM , the effective thermal conductivity n e ff is given 
by 

1 _ 1 

K'eff KO 

where no is the intrinsic lattice thermal conductivity of the material, Rk is the interfa- 
cial thermal resistance and d is the grain size. This expression states that the effective 
resistance is a sum of the intrinsic resistance plus the sum of the resistances for each 
boundary crossed. The relation estimates the reduction of thermal conductivity in a 
material from its intrinsic value due to interfacial resistance. The effective thermal 
conductivity n e ff was evaluated using thermal parameters taken from our atomistic 
simulations. The value of kq cannot be measured directly, therefore, we assigned a 
value of 50 W/ (m ■ K) which was obtained from our previous atomistic computations 


Rk 

~d 


( 7 ) 
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of the bulk lattice thermal conductivity of ZrB 2 [10]. This value is the average of 
the in-plane and normal conductivities obtained from those simulations. For Rk, we 
took the value for a-tilt which has the highest interfacial resistance of our structures, 
3.3 (m 2 • K)/GW. For d we used 6/jm which is typical for these materials. The effec- 
tive thermal conductivity K e ff obtained was 48.6 W/(m ■ K ) which represents a very 
small reduction in thermal conductivity from the intrinsic value. Experimental ther- 
mal conductivity values reported by Zimmerman et al [7] are around 22 W / (rn ■ K) 
at room temperature. This suggests that the boundaries considered in this paper 
have very low interfacial thermal resistance compared to experimental data. This is 
not surprising since our structures are very pristine and current processing methods 
for ceramics typically give much more disordered boundaries with numerous impu- 
rities. On the other hand, this also shows that significant improvement in material 
properties can be obtained by improving the quality of the grain boundaries. 

Next, to consider a more realistic microstructure than BLM, which assumes one 
dimensional conduction and a regular, evenly spaced grain boundary network, we 
evaluated K e ff with FEM computations whose mesh was constructed on an SEM 
micrograph taken from Ref [3]. A FEM thermal analysis requires, in addition to a 
structure, parameters for the thermal conductivity, the specific heat and the density 
of the grains and the boundaries. For the density and the specific heat, we used 
bulk values as measured experimentally. For the thermal conductivity of the grains, 
we again assigned a value 50 W/(m ■ K). Since the MD interfacial resistances are 
significantly lower than what is expected for these materials, a more realistic value for 
<jk was chosen to be 10.0 MW /(m 2 ■ K) [37]. To obtain an interphase conductivity 
kk-, we multiply a K by a representative width and obtain n K = 1 W/(m ■ K) 

The FEM transient analysis was performed as discussed previously. In Fig 10, 
we show the original SEM image, the FEM mesh, and the steady state tempera- 
ture distribution in the material. As we can see, the conduction is not uniform as 
the grain boundary network impedes the flow, resulting in jagged thermal contours 
through the material. The effective thermal conductivity was determined by fitting 
the transient, temperature response at the bottom of the model to an equivalent 
homogeneous, reference material. Remarkably, the result for the effective thermal 
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conductivity of the material is 18 W/{m ■ K) which is very close to the experimental 
value of 22 W/{m ■ K ). Transient temperature traces for the polycrystalline model 
and for the homogeneous, reference material are shown in Fig 11. The trace for a ref- 
erence material with the intrinsic thermal conductivity Ko is also plotted to show the 
reduction in thermal conductivity. In addition to the transient computations, results 
from the UTG and UHF steady state computations are reported in Table II. The 
numbers agree very well with each other and with the transient conductivity indicat- 
ing that this microstructural model is of sufficient size to be considered representative 
of a bulk system. 

It is interesting to compare the FEM result with the BLA4 and also with the 
rule of mixtures using the same parameters. Using the BLM expression and the 
more realistic interfacial parameters, we obtain 10.1 W/{m ■ K) for ft e //. This is a 
reasonable result since the BLM is a series resistance model and therefore should be 
lower than the FEM model which has both series and parallel contributions. Since 
the FEM model is a two phase system we can also apply the rule of mixtures which 
gives for this situation 44 W/{m ■ K). The rule of mixtures describes a system of 
parallel resistances, therefore it should give an upper bound for this model. 


VIII. CONCLUSION 

In this paper, we used a combination of ab initio computations, atomistic simula- 
tions and FEM calculations to study the structure and properties of grain boundaries 
and their impact on lattice thermal conductivity in ZrB 2 . Four CSL and near CSL 
grain boundaries were considered. CSL boundaries are generally low energy struc- 
tures with favorable properties. In particular, we considered tilt and twist boundaries 
relative to the c and a axes. 

Ab initio methods and interatomic potentials were used to evaluate the energetics 
of these boundaries. In particular, ab initio results showed that the c-axis boundaries, 
c-tilt and c-twist, had lower interfacial energies than the o-axis boundaries, c-tilt and 
c- twist. The a-axis boundaries contained a small amount of intrinsic strain relative 
to the c-axis interfaces, but were also much more disordered as indicated by the 
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significant amount of relaxation. Energetics from the interatomic potential did not 
produce as clear an ordering, giving interfacial energies very similar to each other. 
The absolute values of the energies from the potentials were reasonable however. In 
addition, the degree of relaxation of the structures matched very well the ab initio 
results. In general, the new interatomic potential developed for this material seemed 
to give a reasonable description of the boundaries. 

Ab initio computations also produced the electron localization function (ELF) 
which provided information about bonding across the interfaces. From the ELF, c- 
tilt and a-twist were seen to have significant covalent bonding across the interfaces 
whereas c-twist and a-tilt were more ionic. Strong covalent bonding is expected 
to improve the mechanical integrity of the boundary as well as to reduce thermal 
interfacial resistance. 

Next, the interfacial thermal conductance a k was evaluated using nonequilibrium 
molecular dynamics. The Muller-Plathe method was used to produce a heat sink 
and a heat source on opposite sides of the boundaries. This temperature differential 
resulted in a discontinuity at the interfaces proportional to their thermal resistance. 
For our boundaries, c-tilt had the lowest thermal resistance. This result was not sur- 
prising since this boundary had the least reconstruction and is held together by strong 
covalent bonds. The highest interfacial thermal resistances was for a-tilt which was 
both more disordered and more ionic than the other boundaries. In general, however, 
all four boundaries had very low thermal resistances compared to estimates based on 
experimental data. That is not surprising since these boundaries are very crystalline, 
highly ordered, and free of impurities. Modern processing methods typically result 
in ceramics with much more disordered boundaries and with many impurities. These 
calculations show however that low energy structures with very favorable properties 
exist for these materials and might be produced by improved processing methods. 

Finally, the effect of the grain boundary network on the bulk thermal conductivity 
of the material was considered. Grain boundaries are expected to reduce the thermal 
conductivity relative to its intrinsic, single crystal values. Estimates based on the 
Brick Layer Model indicate that boundaries with very low thermal resistance, such as 
the ones we considered, will produce a very small reduction in the intrinsic thermal 
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conductivity. To examine a more realistic situation, we performed a FEM thermal 
analysis using a microstructure obtained from an SEM image and also using more 
realistic thermal parameters. The FEM mesh was constructed directly on top of the 
image. From these computations, a significant reduction of thermal conductivity was 
seen with values much more comparable with experimental results. 
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Ab Initio 

DLB/Pot 1 


7 (rneV/A 2 ) A 2 (A) 

7 (meV/A 2 ) A z (A) o K ( GW/(m 2 • K)) 

c-tilt 

153(369) 

112(238) - 1.77 

c-twist 

157(375) 0.29 

111(258) 0.29 0.58 

a-tilt 

227(1040) 

107(1380) - 0.38/0.31 

a-twist 

212(1230) 

118(1430) - 0.55/0.53 


TABLE I: Energetics and thermal conductance for ZrB -2 grain boundary structures from 
empirical potentials (DLB/Pot 1) and ab initio/ DFT. Units for 7 are nreV/A 2 , A, are A 
and ok are GW/(m 2 • K). Unrelaxed quantities are in parenthesis. 


K eff 

Vertical 

Horizontal 


FEM/Transient 

18 


18 

FEM/UTG 

17.5 

16.2 

16.8 

FEM/UHF 

16.7 

15.9 

16.3 

BLM 



10 

ROM 



44 

Exp. 



22 


TABLE II: Effective thermal conductivity values (W/(m-K)) for the SEM microstructure 
from FEM computations, the Brick Layer Model (BLM), the Rule of Mixtures (ROM) and 
the experimental value [7]. FEM computations are obtained from a transient method as 
well as two steady state approaches using the uniform temperature gradient UTG boundary 
condition and the uniform heat flux UHF boundary condition. Vertical and horizontal give 
the direction of heat transport. The last column gives the average if the result is direction 
dependent. 
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FIG. 1: Atomic structure of £7 symmetric c-axis tilt boundary for ZrB% with Zr (red) 
and B (gray) atoms shown. The boundary contains 5-7 membered ring pairs separated by 
hexagonal rings. This structure is similar to ones proposed for graphene. 
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FIG. 2: Atomic structure of £7 symmetric c-axis twist boundary for ZrB 2 with Zr (red) 
and B (gray) atoms shown. A layer of Zr has been shifted relative to a B layer resulting 
in crumpling. 
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FIG. 3: Atomic structure of 90°, near coincidence, asymmetric a-axis tilt boundary for 
ZrB 2 with Zr (red) and B (gray) atoms shown. Unit cell has two boundaries designated 
Left and Right. A rotated view of this boundary can be seen in Fig. 7. 
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FIG. 4: Atomic structure of 90°, near coincidence, asymmetric a-axis twist boundary for 
Zr £>2 with Zr (red) and B (gray) atoms shown. Unit cell has two boundaries designated 
Left and Right. Significant reconstruction makes this the most disordered boundary. 
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FIG. 5: ELF for the c-tilt boundary. Blue is high ELF and green is low ELF. The 
figure shows a cut through a B plane. High ELF between B atoms indicate strong, highly 
directional covalent bonds form across the interface. 
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FIG. 6: ELF for the c-twist boundary. Blue is high ELF and green is low ELF. The twist 
results in the misalignment of the Zr and B planes. The figure shows a cut through a Zr 
plane where the ELF is low and less directional than Fig. 5. ELF accumulation points 
can be seen forming 6-membered rings around the Zr atoms. 
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FIG. 7: ELF for the a-tilt boundary. Blue is high ELF and green is low ELF. The bound- 
ary is asymmetric. The figure is a cut with maximum ELF in the boundary. Interacting 
pockets of high ELF can be seen due to B pairs. 
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FIG. 8: ELF for the a-tilt boundary. Blue is high ELF and green is low ELF. The bound- 
ary is asymmetric. The figure is a cut with maximum ELF in the boundary. Significant 
and complex B bonding can be seen across the boundary. Formation of 4-member ed rings 
is observed. 
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FIG. 9: Temperature profiles versus simulation cell distance normal to the grain boundaries 


for c-tilt (upper left), c-twist (upper right), a-tilt (lower left), and a- twist (lower right). 
Cell distance is in reduced units. Grain boundaries are positioned at z = .25 and z = .75. 
Temperature discontinuity is related to the interface thermal resistance ax- 
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FIG. 10: Finite element thermal analysis for a ZrB 2 microstructure. The FEM mesh was 
constructed on top of an SEM image of ZrB 2 [3]. The figure shows the original image, 
the FEM mesh, and constant temperature contours from a steady state thermal solution. 
Yellow is high temperature and blue is low temperature. 
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FIG. 11: Transient temperature response at the bottom of the model after a flux has 
been applied to the top as a function of time. The effective temperature of the model is 
determined by fitting the solution to an equivalent, bulk, reference material. A bulk system 
with the intrinsic thermal conductivity kq shows the reduction due to interfacial resistance. 
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